{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4d. Calculate isobars and isopleths\n",
    "\n",
    "This allows you to calculate H<sub>2</sub>O-CO<sub>2</sub> isobars and isopleths for a given melt composition and temperature."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python set-up\n",
    "You need to install VolFe once on your machine, if you haven't yet. Then we need to import a few Python packages (including VolFe). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Install VolFe on your machine. Don't remove the # from this line!\n",
    "# pip install VolFe # Remove the first # in this line if you have not installed VolFe on your machine before.\n",
    "\n",
    "# import python packages\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import VolFe as vf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'1.0.2'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vf.__version__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define inputs\n",
    "\n",
    "The melt composition and temperature can be given in a dataframe, or read from a csv file.\n",
    "\n",
    "In this example it is read from a dataframe,which is from Brounce et al. (2014) with the updated Fe<sup>3+</sup>/Fe<sub>T</sub> from Cottrell et al. (2021).\n",
    "\n",
    "Note the volatile content of the melt is not used in this calculation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the melt composition, fO2 estimate, and T as a dictionary.\n",
    "my_analysis = {'Sample':'Sari15-04-33',\n",
    "           'T_C': 1200., # Temperature in 'C\n",
    "           'SiO2': 47.89, # wt%\n",
    "           'TiO2': 0.75, # wt%\n",
    "           'Al2O3': 16.74, # wt%\n",
    "           'FeOT': 9.43, # wt%\n",
    "           'MnO': 0.18, # wt%\n",
    "           'MgO': 5.92, # wt%\n",
    "           'CaO': 11.58, # wt%\n",
    "           'Na2O': 2.14, # wt%\n",
    "           'K2O': 0.63, # wt%\n",
    "           'P2O5': 0.17, # wt%\n",
    "           'Fe3FeT': 0.177}\n",
    "\n",
    "# Turn the dictionary into a pandas dataframe, setting the index to 0.\n",
    "my_analysis = pd.DataFrame(my_analysis, index=[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll mostly use the default options..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            option\n",
      "type                              \n",
      "COH_species     yes_H2_CO_CH4_melt\n",
      "H2S_m                         True\n",
      "species X                       Ar\n",
      "Hspeciation                   none\n",
      "fO2                       Kress91A\n",
      "...                            ...\n",
      "error                          0.1\n",
      "print status                 False\n",
      "output csv                    True\n",
      "setup                        False\n",
      "high precision               False\n",
      "\n",
      "[78 rows x 1 columns]\n"
     ]
    }
   ],
   "source": [
    "# print default options in VolFe\n",
    "print(vf.default_models)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... but the 'COH_species' must be set to 'H2O-CO2 only' because the isobars and isopleths are calculated assuming the only melt and vapor species are H<sub>2</sub>O and CO<sub>2</sub>O."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# change just the \"COH_species\" option to \"H2O-CO2 only\"\n",
    "my_models = [['COH_species','H2O-CO2 only']]\n",
    "\n",
    "# turn \"my_models\" to dataframe with correct column headers and indexes    \n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run calculations\n",
    "\n",
    "The isobar calculation is run as below. The initial and final pressure, as well as the pressure step, can be specified."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate isobars starting at 1000 bars, ending at 4000 bars at 1000 bar intervals\n",
    "isobars = vf.calc_isobar(my_analysis,models=my_models,initial_P=1000.,final_P=4000.,step_P=1000.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The isopleth calculation is run as below. The highest pressure, as well as the step-size for the H2O mole fraction in the vapor, can be specified."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate isopleths up to 4000 bar, at XH2O = 0.20 stepsizes\n",
    "isopleths = vf.calc_isopleth(my_analysis,models=my_models,final_P=4000.,step_XH2O=0.10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can plot them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'H2O (wt%)')"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAFzCAYAAAAg407BAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAYsVJREFUeJztnQWYE1fXxy/u7u7ubsWKbHEoFF6KuxSKFW9xSoEixd1KKS2uxd2d4u66UByKzvf8z/vefJPsbjbZnclIzu95ht1NQnJnktxz75H/iaAoiiIYhmEYJgQihnQHwzAMwwA2FAzDMIxb2FAwDMMwbmFDwTAMw7iFDQXDMAzjFjYUDMMwjFvYUDAMwzBuYUPBMAzDuCWy+7sZ8OnTJ3H37l0RJ04cESFCBKOHwzAME25Qa/3ixQuRMmVKETGi+z0DGwoPgJFIkyaN0cNgGIbRnFu3bonUqVO7fQwbCg/ATkJe0Lhx4xo9HIZhmHDz/PlzWgDL+c0dbCg8QLqbYCTYUDAMYyc8cadzMJthGIZxCxsKhmEYxi1sKBiGYRi3sKFgGIZh3MKGgmEYhnELGwqGYRjGLWwoGIZhGLewoWAYhmHcwoaCYRiGcQtXZvtAeOvdu3eO4+3btyH+LX//8OGDiBcvnkiSJAkdiRIlEpEj81vFMIwx8OyjI927dxfjxo3TpMQ+QYIEDsMR0pE5c2aRPn16VrhlGEZT2FDoSHC7AEzi0aJFoyNq1Kh0uP4eKVIk8fTpUxEYGCj++ecf2pXgJ44LFy64fU3sRPLnz+905MyZk56bYRgmLERQMAsxoaosYgJ+9uyZV6KAeDzcSa5GwBvghoKBgNFwdzx48EBcunRJvH//PshzRIkSReTKlcvJeOTLl0/Ejx/fq7EwDOOf8xobCh0Nha9BfOPcuXPixIkTTgd2J8EBN1XRokVF1apVRZUqVUTSpEl9PmaGYYyBDYWfGorgwNt748aNIMYDt7m6xIoUKSKqV68uqlWrJgoUKMCxDoaxMc/ZUGiLlQ1FSMCdBYOxY8cOsXbtWnH8+HGn+1OkSEE7DRiOihUritixYxs2VoZhtIcNhcbY0VAE1+51/fr1ZDS2bNkiXr165bgP8ZWyZcvSTgOGI1OmTIaOlWGY8MOGQmP8wVCoQQB+586dYt26dWQ4rl696nR/9uzZRePGjUXLli1p58EwjPVgQ6Ex/mYo1ODjgZRcaTT27NlDmVgAGVw1a9YUbdu2FZUrVxYRI3KhP8NYBTYUGuPPhsIVXIOVK1eKmTNnir179zpuT5cunWjTpo1o0aKFSJkypaFjZBgmdNhQaAwbiuA5ffo0GYwFCxY4UnCxy6hRo4Zjl+Ft3QjDML6BDYXGsKFwz5s3b8TSpUvF9OnTg+wyWrduTbEM3mUwjLlgQ6ExbCg858yZM7TLmD9/vtMuA9lS33zzDaXacn0Gw1hrXuPoI6MpkAoZP348pdvCJfXZZ5+Jjx8/ilWrVpErqmTJkmLTpk0UJGcYxhqwoWB0IUaMGKJJkyZi9+7dtMvo3LmziB49ujhw4IAICAgQpUuXpnoNNhgMY37YUDC6A/XaCRMmUD1Gly5dSBwRsYxKlSpRId/27duNHiLDMG5gQ8H4DBTnwS0Fg4EdBgwGdhyff/65KFeuHBX5MQxjPthQMD4HGVDYYVy5coUC3JAIgZGAsYDRgPFgGMY8sKFgDCNVqlRi0qRJ4vLly6J9+/bUNwNuqDJlypBbat++fUYPkWEYNhSMGUiTJo2YOnUqNV5CoR46AyLQXapUKQp8uyrbMgzjW9hQMKYBBXoo2oPBQKEe6i+QSlu4cGHx7bffUr43wzC+hw0FYzrQeQ9FexcvXhQNGjQQnz59EhMnTiTV2kWLFnFKLcP4GDYUjGnJmDGjWLx4sdi8ebPImjWruH//vmjUqBFVd58/f97o4TGM38CGgjE9MAx///23GDZsGBXtbdu2TeTNm1f069dPvH792ujhMYztYUPBWALUXPTv31+cPXuWOu29f/9ejBgxgor5Vq9ebfTwGMbWsKFgLEWGDBnEmjVrqCdG2rRpxY0bN0StWrWogdL169eNHh7D2BI2FIzlgPosjAN2F3369KF0WhgP7C5+/PFHauXKMIx2sKFgLEusWLHI/XTy5Emq6kZfDLin8uXLJ7Zu3Wr08BjGNrChYCwPdhIIcC9cuFAkS5aMenwjAN6xY0cOdjOM1Q0FVoNFihQRceLEEUmTJhW1a9emL7maf//9l/SAEiVKJGLHji3q1q0rHjx44PSYmzdvUoAzZsyY9Dw9e/YUHz58cHrMjh07RMGCBSkomjlzZjFv3jyfnCPjO3cUUmeRNgsDAVDtjWK9EydOGD08hrE0hhoKCMHBCKBHAXLlkcmC5javXr1yPKZbt27kf16yZAk9Hg1xvvzyS8f9aIoDI/Hu3TvSBkJnNRiBAQMGOB5z7do1ekz58uVp0ujatStV/m7cuNHn58zoS/z48cXkyZPpvU2ePLk4d+6cKFasmBg7diwV7jEMEwYUE/Hw4UOU3Co7d+6kv58+fapEiRJFWbJkieMx586do8fs37+f/l6/fr0SMWJE5f79+47HTJ06VYkbN67y9u1b+rtXr15Krly5nF6rQYMGSkBAgEfjevbsGb0mfjLWITAwUKlZsya9dzgqV66s3L171+hhMYwp8GZeM1WMQmr5JEyYkH4ePXqUdhnwN0sg44C0yP3799Pf+JknTx7yTUsgJId+sOisJh+jfg75GPkcriBrBv9ffTDWI3HixJRGCxcUOu5BNwqfFa67YBjvMI2hgFsALiEohubOnZtug2QDehXAnaAGRgH3yceojYS8X97n7jEwAMiUCS52gqbj8oC6KWPd2AUkzLHoyJ8/v3j8+DGl1nbo0IED3QxjNUOBWMXp06dJ28do+vbtS7sbedy6dcvoITHhJEeOHBQL69GjB/09bdo0DnQzjJUMRadOncTatWupaU3q1KkdtyMYiSD106dPnR6PrCfcJx/jmgUl/w7tMXHjxiWXhCvIjMJ96oOxPnhff/75Z3JBoS0rB7oZxgKGAnLRMBIrVqygPHjIM6gpVKgQdT1TF08hfRbpsCVKlKC/8fPUqVPi4cOHjscggwqTO/Lr5WNcC7DwGPkcjH+B7nkQGYQLCgsR7DK++OILyqhjGCYYFAPp0KGDEi9ePGXHjh3KvXv3HMfr168dj2nfvr2SNm1aZdu2bcqRI0eUEiVK0CH58OGDkjt3bspoOXHihLJhwwYlSZIkSt++fR2PuXr1qhIzZkylZ8+elDU1efJkJVKkSPRYT+CsJ3vy6dMnZdq0aUqMGDHo/U2UKJGyatUqo4fFMD7Bm3nNUEMh0xZdj7lz5zoe8+bNG6Vjx45KggQJaLKvU6cOGRM1169fV6pUqUJf+MSJEys9evRQ3r9/7/SY7du3K/nz51eiRo2qZMyY0ek1QoMNhb05e/YsfTbk569///7Kx48fjR4Ww+iKN/NaBPwT3E6D+X+QHYXsJwS2OV5hT5ASjSSGcePG0d9wS/3666+kGsAw/j6vmSKYzTBmCHQjqL1gwQL6fdWqVRTDunLlitFDYxjDYUPBMCqaNGkidu3aRVlRKNiEFhkr0TL+DhsKhnGhaNGi4siRI/TzyZMnVMU/YcIEytJjGH+EDQXDBEPKlClJhBI7DAhPdunSRbRp04abIjF+CRsKhgmB6NGjkxoxivQiRowoZs+eLT7//PMgxZsMY3fYUDBMKFpRKMhbv349ZYhAyh7SH9COYhh/gQ0Fw3gA4hSHDh0S2bJlE7dv3xafffaZKXTJGMYXsKFgGA/JmjWrOHjwoKhSpQp1XmzYsKHo168f60QxtocNBcN4AdxP6LjYq1cvhyQ9ivO4ZwljZ9hQMIyXRIoUSYwcOZIqt1GcB+VjFOexHD1jV9hQMEwYady4sdi9ezel0p49e5aabp0/f97oYTGM5rChYJgw8ueff9JuYt68eRTkxo6idOnSVKzHMHaCDQXDhJE5c+aIIUOGUD8U7CzQP+XRo0eifPny1F+FYewCGwqGCSNNmzYVDRo0EFWrVhVJkiShDo0oyHv58iVlRi1fvtzoITKMJrDMuAewzDjjKUibbdSoERkJVHNPnz5dtG7d2uhhMUwQWGacYQyU/UDsAsYB9RXQhxo1apTRw2KYcMGGgmHCwLp160hZNqT02RkzZog+ffrQ37179xY9e/Zk9VnGsrChYBgvuXHjhqhevbpInjy5ePXqVYgaUSjGg6AgwM9WrVqJDx8++Hi0DBN+2FAwjJfcu3dP5MyZk/pVxIoVy+1jISg4d+5c2mXgZ7169SiOwTBWgoPZHsDBbCY4Xr9+LWLGjOnRY9FaFRlS6GdRtmxZ+hufKYYxCg5mM4wP8NRIAOhBbdy4UcSJE4caIqHW4uHDh7qOj2G0gg0Fw3gBdgRh3YRjJwEjkTRpUnH8+HGSKr9+/brmY2QYrWFDwTBegErsdOnSUbe7sFCgQAGxZ88eeo5Lly6Rsbh8+bLm42QYLWFDwTBeAGkOaDpFiRIlzM+RJUsWsXfvXgqI37lzh6q5r127puk4GUZL2FAwjBds3bqVaiiqVasWrudJlSoVGZ3s2bOT4YGxuHnzpmbjZBgtYUPBMF4GsKHtlChRonA/V7JkychYYIeBWAWMBXYYDGM22FAwjIGkSJGCjEXGjBnFlStXyFjcv3/f6GExjBNsKBjGA1BRjWK5iRMnal4wlzp1ajIWadOmFRcvXhQVKlTg1FnGVLChYBgPOHDggFi2bJkYNGhQuALZIYEsKMiUI3aBbnkVK1YUjx8/1vx1GCYsRA7T/2IYPyNDhgykAvv+/XuS49ADuJ+ws0C9BZohVapUiYLnCRIk0OX1GMZTWMLDA1jCg/El586dI2MRGBgoihQpIjZv3sxyH4zmsIQHw1iYHDly0E4CmVWHDx+mLKsXL14YPSzGj2FDwTChgMkaLqF379757DXz5MkjtmzZQm6nffv2kax5SJLmDKM3bCgYJhRGjx5NmUjoL+FL8ufPLzZt2kRugV27domaNWuKN2/e+HQMDAPYUDBMKKBBEYT8AgICfP7ahQsXJtXZ2LFj066mdu3a3M+C8TkczPYADmYz6H8NIkY0Zm0FIUEYKvTAgHzI8uXLRdSoUQ0ZC2MPOJjNMBoDA2GUkQBQmYXGVIwYMegnmiAhVZdhfAEbCoZxw6NHj4RZKFeuHHXGixYtmli5cqXo0KFDmHtjMIw3sKFgmBC4ffs2xSaKFStmmtU7ivD+/PNP2t2gJ8bgwYONHhLjB7ChYBg3sh1YsaMSWw/ZjrCC7KcpU6bQ7zAUs2bNMnpIjM1hCQ8dQTMarEoRCP348WOQn663AeTNJ0mSRCROnJiO6NGjG30afgtEAO/du2dKgb527dpRH4vhw4eL9u3bi5QpU1JhHsPoAWc96Zj11L17dzFu3LhwvTbSIqXhUP/EkTVrVpE7d27SCNJLf4gxL/jqtmjRQsyfP5/6ZOzYsYMkPxhG63mNdxQ6Av82mtJgEodPGT/d/Y5dxZMnT0jjB0FUSFu/fPmSDnetMrHrgOxDrly5yHDIn5CtNjJTh9GXCBEiiJkzZ9KuB4V5SJvdv3+/yJQpk9FDY2wG7yhMWkeBtwWvB4OBQxoP+ROTA8TjIEkdUrUudiPoywyjUbx4caouxu6DCR3IiV+9elV07tzZ9Kt06EBBRPD48eMic+bMJPmBHSfDaDWvsaGweMEd4htoo3n69Glx5swZ+onj/PnzwWbqpE+fnrqowWiUL1+eOqwxzuArgQkXhgJpqLVq1RJmBwuHEiVKiBs3boiiRYtSFXesWLGMHhZjYthQ+JGhCAkYicuXL5PxOHnyJPmvkcUDd5Ya7DhgNGA8kKcfP3584e/gK7Fz504qbBs4cCDtzKwAFgelSpUS//zzD4kIrlixQkSOzN5lJnjYUGiMFQ1FcCDWASkISFhjxQlXhfrtRzyjYMGClD3zn//8h+IejLWA2wmGH3pQbdu2FdOmTaNYBsO4woZCY+xiKFxBq02snKXhwIrUVb3066+/JqORJk0aw8bJeAd2EnXr1qVFwNChQ8X3339v9JAYE8KGQmPsaihcuXv3LvVAWLJkidiwYYOTmwpaQzAaqC2wc6AU7zVW4cggQvaYVZk8ebLo1KkT/T5nzhxKo2WYMM9rMBSMe549ewZjSj/9hUePHinTp09XypUrp0SIEIHOH0ekSJGUKlWqKAsWLFCeP3+u2I0lS5bQeWbLlk2xOr1793a8Z3/99ZfRw2EsPK9xkj0TLGjDCR/39u3bqQJ4zJgxolChQpRl9ddff4mmTZtSnQjcUoh72GVjisr4L774glw3VufHH38UjRo1ovcMO8GjR48aPSTGqigGsnPnTqV69epKihQpyLKtWLHC6f5mzZo5VrLyCAgIcHrM48ePla+//lqJEyeOEi9ePKVly5bKixcvnB5z8uRJ5bPPPlOiRYumpE6dWhk5cqRX4/THHUVIXLhwQRk0aBCtuNXvS+HChZWFCxcqb9++NXqIjAq8HxUqVKD3KFmyZMrVq1eNHhJjEiyzo0AP4Hz58pE/NSSwukOOuDx+//13p/uxYkIK6ObNm8XatWupZSRWwmo/XOXKlUW6dOloRYW2liimmjFjhq7nZlcgG4KUURT7HTlyRLRu3Zpkr/F748aNRYYMGWglayZ5bn8GzY3Q5AjfswcPHtD3id8bxmsUkxDSjqJWrVoh/p+zZ8/S/zt8+LDjNvhi4VO/c+cO/T1lyhQlQYIETitd+G698UHzjsI9Dx8+VIYNG+bYGeKIHj260qZNG+XMmTOKVTh//rzy6tUrxY7g+5A2bVp6bxB3evfundFDYgzGMjsKT0ChGHzh2bJlo0YtSOmUQNcGBWLoKyypWLEi1QMcPHjQ8ZgyZco4tY1ES8kLFy6QrlJwvH37lnYi6oMJGWRB9e/fnyrEf/31V6rFQB4/dIiQOYTrjSwqqZBrVuDHR2wGcRm7AXVZFBCieBDfKQhWMoynmNpQYJu8YMECyvMfOXIk5fxXqVKFgnPg/v37ZETUoBI1YcKEdJ98TLJkyZweI/+Wj3FlxIgRlDYmD64h8AwYY7if4IaCC/DLL78kow3BOrxvMBowJPL9MxNIEYRmEira4aaxI9D8WrhwIf0+adIk7mPB2MNQIKMGTVry5MkjateuTTGIw4cP04pIT/r27UsThzyQ9cN4DiqBS5cuLZYtW0YyIt26dRNx4sShgj5kS6GQDy09zZQphQUBFHqx08RCw65At2rIkCH0e8eOHamSm2EsbShcgfIp+jFg8gHJkycP0lQGRWLQusF98jEI4qmRf8vHuILgLApQ1AcTNhDcHjt2LDVwwk4NrkKIFsLwlyxZUnej762B8weJbrgJkf6L3RN2fXhvGMY2hgIfaMQopOIp1DKfPn3qlB8OKQr4wtHnWD4GbhC1kioypBDzQM484xtgbPv06UOKrNixodEORAqhYIsYhpE5/mba2fgCuAPnzZtHO3UsmmC0Q5KqZxgiPFHzf//9Nzz/neodjh8/TgeGMnbsWPr9xo0bdN93332n7N+/X7l27ZqyZcsWpWDBgkqWLFmcXveLL75QChQooBw8eFDZs2cP3d+wYUPH/U+fPqX88SZNmiinT59WFi9erMSMGZOqjj2Fs5605969e8o333yjRIkSxZEpVa9ePeXcuXOGVGNnz56dPn/+BGoqEiVKRNe+UaNGyqdPn4weEuNDvJnXvDIU69evV5o2bapkyJBBiRw5shIxYkQqdCtTpgylR8qUVE/Zvn17kII6HEiLff36tVK5cmUlSZIkNJmkS5eO0i3v378fpOAOhiF27NhK3LhxlRYtWrgtuEuVKpXy008/eTVONhT6ceXKFaVx48YOmRB8plA0icWCr5CFnd27d1f8jW3btpHEB85/9OjRRg+H8SHezGseiQJCjbJ3796UFQIJajRGQbpdjBgxKB4An/Pu3bspFbV58+akWGkn4Th/EQU0klOnTpHK6erVqx0ZVOguh+I+BML1BO7LjRs3OtrI+hvIgMK1hksKKbTINmTsz3OtRQGLFy+urF27Vvn48aPbx92+fZuK2ey2hecdhe/Yt2+fUrZsWcfuEjvApUuXsltER3BtW7VqRdcbMjiQaWHszzOtdxT+Du8ofAs+khAexCoXwW+AnSxWvsiiYrQHRabocoh02ezZs1OiAT7zjH3xZl7zOutJfnEZRs80VRgGuDThjooSJYpYv349tW1Fiu27d+80ey30bJg/f754/fq18GeQEo66F9SQoN4FvUfMWBjJGIS32xUEHdOkSUMByFmzZimXLl1S7A67nowFmVDly5d3uKNy5Mih7NixI9zPCxcLng/JEsiOYxRSc5bXuV+/fkYPh7Gq1hOqlLGqQyB71KhRpCaaOnVqUnFlSQBGD+AKgYwL5D+QJAHl2nLlylHiRGBgYJifF7pHAwYMILVhdrP8FyQToN4FQAX4jz/+MHpIjBkIr1W6ePEipRfKdFk7wjsK8/DPP/8o7du3d6TTQhl45syZoSZaMN7Rs2dPur4xYsRQjh07ZvRwGKvtKODLhchbv379SIIhb9684uTJk+Trhe49w+gJqumnTp1KQVeI90EBuE2bNtTTG751JmwcOnSIAtoSeA2QJouKbVRuu0rlMH6Gt1YI/tykSZMq3bp1U1atWkUrPLvDOwpz8v79e2XcuHFUbClXv5MnT/YolRY9TFDNzzuR//YTwTVEUev169cdt+O7DaUDXNvSpUtz90KboeuOAtkoyIZYvHgxHUuWLBEXL17Ux4oxjBsgKd+1a1eKWaCLIVa/33zzjahevXoQIUhXoKAKPbCff/5Z+DuXLl2i9Ej04lBL6mP3BpVfFDyioBYqwIx/4rWhWLlyJbVSRCMaCO7BDQVJ6VSpUlFAm2F8DZIpUHfxyy+/UJonUmkheLdmzZoQ6zSQBopgNlchC3IhQ5H5zz//pOpsNTly5BCLFi2i36dMmcLBbT8lzAV3+G/Hjx+nbmA4IIGA2yDzbTe44M46oPYCC5a///6b/m7Xrp0YM2aMiBUrVpDHwicPqRDUbTChS5MjCwq7Cyj9ZsmSxeghMWYuuENvATQTwjYVW/fff/+dUmRRrBOeVEWG0QLoNSEw26NHD/p7+vTp1JoVXfdcwe7Dn40EXMaQ4PeEwYMHU0th6L199dVX1OqW8SO8DYAULlxY6dGjh7JmzRq/KVLiYLY1gTQ9tKLw3iF9e/jw4RSQ9YcEDE+oU6cOXRtcF0+AOjTUnPF/kKLMWBvdZMb9FTYU1gUy9PXr13dUG+fJk4dktevWrav4Mx8+fKB+IJDeP3v2rMf/b+PGjY4alkWLFuk6RsbiooDYdsIHjNxqdJNTA7eU3eAYhbXBRxxV3aj1gesElCpVSuzZs0f4O3AXe9sSAPpbw4cPp2QAxCvgembsPa95bSiQ7dSkSRNqSRrkySJEsKWQGBsKe3Dt2jX67O7du5f+7tixoxg3bhwFtBnPQcJKxYoVxc6dO6noEX1oIOnDWAtdg9mQfq5fv764d+8e7SbUhx2NBGMfIFGOyW3QoEG0qEG6J6S18Vn2p0ke+lbhOWfUryBlFjsRqDJwfYUf4K1fC61PL1++rPgTHKOwH2jEhSY9eF9TpEih7N27V/EHZsyYQeecNm1aqmwPD5s2beJ4hYXRtTK7Xr16YseOHfpYLYbRERTXoXIbq+lq1aqJw4cPU+tT/A01Wuww7N7HC+nDxYsXp4p27AzCQ6VKlSheAaDAe+HCBY1GyZgNr2MUEAVEHjW2nah+RVMZNd9++62wGxyjsIdMBYKumByRhAF5CvDy5UvRqlUrqkoGLVq0IIMRPXp0YVfwlYebOLyGAuB5EK/A4hECoeiMx/EKa6BrMHv27Nmiffv29EVC0Z26YAm/27EDHhsK64MqbPS0gMJs9+7dne7DVwDV271796ZYW+HChamANG3atIaN10pgR5Y/f34ywFDynTFjhtFDYrSe17z1ayVLlowKdPxJdZNjFP7B5s2blUSJEtF7nThxYmXbtm2KXRgyZIgybdq0cMcl3F07Ga9YuHChLq/BWKjgDo1iOJjN2JVr164pBQoUoPcbhXljxozxSLbczFy9epXaA+CctGghGxIDBgyg14gVK5Zy/vx53V6HsUAwu1mzZqwgyViuveekSZPE3bt3Q31s+vTpqc6iadOm5H+HZtTXX38tXr16JawKlJ0hp47zKFu2rG6vg7RbJAXgWiGOCdl3xh54HaNAsHrBggVUaIPglWswG6KBdoNjFNamQoUKYtu2bTRZSrHA0MDXYvLkyVQjgNoDCAuuW7dOJE+eXPfxWj1eUaBAAeoH0rp1azFz5kyjh8QYUXB36tQp+iBAtx6SzpAal8eJEye8fTqG0R3IyqDnQp06dTz+P0jMgOQHDAwy/I4dO0b9V6zUbvV/rmWfvmaKFCnEb7/9Rtdv1qxZYuHChT59fUYnNHJ32RqOUfg3ly5dUjJnzkyfAcTodu/erViB5cuXK2XLlqW2r75m4MCBjnjFuXPnfP76jMExCjW3bt2ig2HsTObMmcW+ffuoUO3JkydUN7B06VJhZrCTgFQJJEvQztTX/PDDD6J8+fIcr7AJXhsK+GvxIYBvC4E/HPgdFZrv37/XZ5QMEwagFLt582ZNPpdwP6EOo1atWlSTAb0zCAqaFbh+0Aq2Q4cOolevXj5//UiRIpEeVLJkychFjRoVxsIoXoKGJUmTJqWc7JMnT9KB35MnT27bZibserImyOfH+1ayZEnN+zjI/hZdu3b1q5oib9mwYYPjWqGXBeMndRRx48ZV1q9fH+T2devW0X12hA2FNcECBh3Zvv/+e02fF3UVo0aNckyAaIL0+vVrxSyY7XMqDSvEFx89emT0cBhfxCjQZxjupuAknFnXnzET7dq1o9oJrd0ecOv07NmT+sXjMw+5DwjkBdejxYieG6lTpxbfffedaVzBo0aNEtmyZaPUWbjC7C68aEe8NhRIGRw6dCj5aSX4HR2vcB/DmAkI36ETmx785z//EZs2bRLx48enIj10zcNEbSSLFy+m2Az6RGgh+qcFMWPGpDRZjGfJkiWUPsvYvOAOuegI6mFngaI7gA/lu3fvqLBJzfLly4Ud4II76/HPP/+IhAkT+uS1zpw5I6pUqUIZgAjerl27loQFjQBf540bN9KuApLiZmLYsGGUCIPvEFopp0uXzugh+TXP9VSPhQyzp8ydO1fYATYU1gLS4UmTJiUZ/PXr15PKsd7AxVW1alVaNGEFvXLlSnJHMc4Zk2XKlKHWqZASwYIT2VGMDdVj/REOZluLrVu3KhEjRlQyZcrkU0E/fD4qVapEn5Vo0aJRFz1fcebMGeXNmzeK2YGgKIrwcI1Gjx5t9HD8mme+KLiD9vzu3bvpwO8MYxbQBxsrfCkl4SuwKkPtQu3atSluBzftihUrdH9dFLOhe1+OHDnIDWZmMmXKJMaPH0+/9+/fn1xQjPmJGJbtSpMmTUiREttHHPi9cePGtIVhGDOAWEGxYsV8/rqI3aFbXoMGDSjrCFXJCDDr3b0PDZegdovsQ7ODjoLQ30Jcs1GjRuLff/81ekiM1oYCHawOHjxIAbunT5/Sgd+PHDlC6YgM4+9AURm7GSlVjslw/vz5ur0eVJzRrxrfQ8RHzA52eVCVRRwJVduy7zZjYrz1a8WMGTNYUbRdu3bRfXaEYxTWoV27dkqTJk2UU6dOGT0Uqthu06YNfXbQ/W369OlGD8lUrF692nFt7NRN0CroGqNABgki5a7gNtmwnmGMAAJ0v/76Kx2vX782ejgkxT99+nTRuXNnSlvFjnvixImaPT8yrLZv3y6sSo0aNchDgWuDhmjwTjAmxVsrhFVRxYoVlXv37jluw++VK1cmyQQ7wjsKa4AV/N69e0myw0ztSzGW7777ziH5AfkPLZ7zs88+o+f75ZdfFKvy4sULh4R7o0aNjB6OX/HMi3nN6zoKNC26fPkyZXWkTZuWbrt58yYF8bJkyeL0WDR7sQNcR8GEF3zNBg4cSKoGYMiQIVR8FlYQAIaMCLpNItMJBXZW5cCBA1TVjoA8Av9IBGAsXnA3ePBgjx+LL4YdYEPBaAWkbmTwFumhMBzhSeGFuwYSIlYH/bZxLXAuCHAjk5KxsKHwR9hQmB/0s0Yr3oYNG4qMGTMKMzNmzBgS7QPdu3enXt6+rPcwI0glRrtaZE+iMRRkSBDjYSzUM5ttCWN2pkyZQit1K/Ro7tGjh5g0aRL9PnbsWEew2xOgJ4Wg+IMHD4TdUorx3sWIEUNs2bJF06A/owGeBD1y5Mih/P7778rbt2/dPu7ixYvUvGjEiBGKneBgtvn59ddflYCAAOX8+fOKVZg5cyalhsoGSJ4E4Bs0aECPr1q1qmJHJk+e7JBAgSwJY6FgNsS7oOl/9epVEjqDMmbKlClF9OjRqYfw2bNnxZ49eyioBqnxfv36BZtCa1XY9cToxbx58xxCm/DThxYDROC3a9euYtq0aSJ//vzCbmA6grjihg0b6PxQ3Mt9biwWo4Ax+OOPP0jf6caNG6QxkzhxYsqECggIoApUO9ZSsKFg9ARuKLifwOjRox3xi5DAV9bOMQ00OILyLxpB9e3bV/z4449GD8mWcDBbY9hQmJfAwECxY8cOUa1aNUvIV4TEiBEjaCcOsFtwlcNB6qg/BXfRy6Zu3bp0zthVGNXfw8481zqYzTBmBQJ89evXJ0NhZbByxgHQLlQdlMfKOmfOnGLWrFlkMPyBL7/8kjoI4nwhIggBQcY4DDUUu3btojJ+xDuwlUazFzXY7MBvmyJFCsqGQNoclDJdO5nB5QWLiBxsfKjQuEYNpIxLly5NMZU0adJQD19fgAwVpPu5HkePHnU6MD481gyyE1YD/mv0cK9Vq5awOrKdMD73zZs3d3wfJkyYQKJ/yATyF0MhzxuSQfh+jBw50ujh+DeKgaxfv17p37+/snz5coq+r1ixwun+n376SYkXL56ycuVK5eTJk0rNmjWVDBkyODVo+eKLL5R8+fIpBw4cILFCyAE0bNjQcT8i+smSJSN5gNOnT1P2VowYMbwSaAtr1lOPHj0csg2eHtGjR1dSpkyp5MmTRylbtqzy5ZdfKq1bt1Z++OEHZcGCBcr+/fuVR48eeTUOu4NsodAy8qwkQ9KsWTP6LESNGlXZtGmT8u7dO2XcuHEkvOlv/Pbbb3QtokSJQt9fRju8mddM0+HO1VDgy588eXKnLlhPnz6ltDlM9uDs2bP0/w4fPux4zF9//UUph3fu3KG/p0yZoiRIkMBpIundu7eSLVs23Q3F0KFDlTRp0jgdqVOnDnIkSZJEiRw5slcGBedUpEgR5euvv1YGDhyoLFy4UDlx4oTy/v17r8bImA+8h3Xr1qX3GYua4NSa/QXMA9WrV6drUaxYMeXDhw9GD8k26Kr1pBdwPaEbGLqDAaTiohvW8ePHndIA0SgJf//yyy9izpw5VLyEFF11X164mJYsWUIdxtATAEEbtVsLipvogga3VXBZWtCxwiHB/4fLSs9gNt6GFy9e0Jjgk8ah/v327dvkdsNx586dEJ8nVqxYFPhD0x552FEOAdcL6di5cuWyZQYQfPKVK1cWO3fupM/ctm3bRKFChYQ/gs87YjT4HqJAsVu3bkYPye+C2ZGFSbl//76jU5ka/C3vw080P1ETOXJkkTBhQqfHuHb9ks+J+4IzFMhA8UbTSgsw2eHNwgGfuzsQy7hy5YrDcOC4ePEiSVjA2GBywSGBoShevDgZDRhaGBKrZ9BAYhtp2QULFhSHDh0SkSJFEnYC7/G5c+ccX2SknyOmhwnT38DnFzInbdu2JX0sdMfDIpLxHRG90WLp1auXyJw5syhatCit5tVAUsAuX1Zkn+DLKQ8Ems0E0kCRZ47MEBRCIhsGkwh2VhBUmz17Nn2p8uXLRwYBK7Jly5bR+wdjAUOJ1rWLFi0Sjx49ElYEAU4oFkPXyS6fOzUwfjD6WAhhJ4FdJZI5sEDwR1q3bi3Kly9PtVuyhwXjQzz1Z8EPjqAwYgYIQCPI3LZtW8f99+/fp9iAVjGKK1eu0G3Hjx93elyZMmWUb7/9ln6fPXu2Ej9+/CD+3UiRIlGAHKDbWa1atZweg25aeO5//vnH9hIeL1++VHbu3KmMHDlSqV27thI3blynWAfeM/h+Bw8erBw6dIiCqVYBMavr168rduXatWvKkSNHKHkhd+7c9H6lT59euXXrluKPXL58mWI2uA7cLdCkwWxkE61Zs8bx96VLl+i25s2bU8AJhiJixIiaB7N//vlnx204oeCC2fgySTZu3BhsMBuZI5K+ffv6JJhtRnAdYDgQ0M+bN2+QIDkC68iy2rFjh6WMht1BczDZ4Cd79ux+m/k2duxYugZx4sTxW4NpakMBS44Vjprbt28rWbNmpdRTTMzeGgp0t8KOAQcGjA8Bfr9x44YjPRY7hlWrVil///037QyCS48tUKCAcvDgQWXPnj1KlixZnNJjserETgg7C6TXLV68mHp7+yI91grgywZxOqTh4sunNhrI0oJBMUP/aTX//vuvYlc2bNhA4prBgd0T3hO8N+hup/4e+AvIekK2H65BtWrVTNXJ0GroYigwQW/ZsiXI7TAQMBaVKlXy2lBs37492NRP5JEDfAhQP4CJHjuJChUqKBcuXHB6jsePH5NhiB07NrlVWrRoQQZIDWow8MXCc6RKlYoMkDfY2VCoQQox3uNWrVoFcVFh94EWnkav4l6/fq0kTpyY3Gh47+0EdgmJEiWi+gm0dA0OLHbg9sV7Ur9+fb/c9eEaYK7BNUCdBWMiQ4HJo2XLlsHeh50FtsXhcT2ZGX8xFGqwWl2yZAlNyCh2Usc0ypcvT7LeRqzs161bR+NImzat7SbJmzdvUu/5XLlyOblKXdm6davjPenZs6fijyCmhvOHYX3w4IHRw7EkuhgKbHuxLQ4J7CzmzZun2BF/NBRqsHKfNm2aUrp0aaddBnZ6SHK4e/euT8cDNyQKK+0IdtGBgYGhPg6GWr4PkyZNUuzOkydPKCFDLg5gSGWMDT06GD+pzDYz/m4o1CBONWTIEJIZkRMVVreoEIeMCuM7hg0bRtcfO3nE8ewclyhRogSdK2RxJEhiQYYjbofMD2MiQ/Hnn38qderUoe0xDvwOF4WdYUMRFKzokBhQsmRJp11G0aJFSU7ELtpLvmDq1KnK8OHDvb5m2H0gQ01KfSC92a5gB4WMPMjUqEGyBc4/RYoUtOtgDDYU2PIheAYfNVJLkYGEA4FsrGiw/bNrBgIbCvdgZYcEBARh1W4prHi1vGY1atRQevXqZSufNNLKkYiBa4bJMCwGG5l/+P9JkyZVrl69qtiV58+fB5vcgDkI5484KmOwoUDqasKECZ1qKSTY9uI+KFzaETYUnoEJHEKIarcUali0MBjIdMHzQTzRToYCiyvswLDoCmtwHhNo/vz56fpgEWeHbLBXr17RosA1gzE4IJooP2+bN2/2yfjsgC6GArLXqIQOiVmzZtFj7AgbCu9XuVgdY9LSymDALYNq+x9//FHz8doBJJPYqcYCLm1ZK+EJ33zzjaNy3RPjwij6GAr0SZCFcCFlReExdoQNRdiDkIsWLaJKYrXBwK7D368lXCbuUmDDgrrGAq5gK6cPIzECEvwh1ZMEt6tCyjTOXUr8MAYYCnzBUbjmLmXRVXfJLrChCB9sMILStWtX0m+CooCWqGss4LqxMt7W6UC+R9b6QKWBMcBQVK1aVWnfvn2I97dr106pUqWKYkfYUOhnMFBljawfdw2XUFwF12ZwwUwrgvNAsB/n7642KayoaywmT56sWEW8smnTpuEWeYT2nNTDsrr7zZKGAltArFS++uorWgXhyaGjhNac9erVo/vsasXZUGhvMCDsqI5hYHUdXCASCr8ym8pVSdjKoKhOTwVUdY3F6tWrFbMD6R2MF22Nw+Myw+cFYqJ4Lsj/MAbUUSCYiBUgPnzqA2X0S5cuVewKGwp9gI9+4sSJlDEnDQZSYNV6XliMoCIXixEmbDUWEMF05zY2A5ABKlSoEC08wwvqunDeWGC4asMxPmqFis5bGzdupK5qIGvWrNSyEc107Io3LQMZ70HL1yFDhojJkydTK1t0KezcubP44Ycfgu1AaFVu3LhBTaRKlizpk9dDs7Fq1aqJzZs3U9fEw4cPi8SJEwuzgqlIi7a2eJ6qVauKDRs2ULOnTZs22bJdrk/ntVBNCcM7Ch9x7tw5SoeUuwvsVOFjdxe/sNIKPyAggAKtvowboKYiU6ZMdD0h5qh1plV44jSIaaIGQq8mR1CLxnnL/jWMD1xPyKbIkSNHsE8K90DOnDmVXbt2KXaEDYVvQYAXnydpMCAVo4VLwkgQWEWgFSnk58+f9+lrI21WVn937txZMQNQvZU9T/RSIZYKs4hZYI5ifGAo4DtGdXZI/PLLLyRJbUfYUBhTmSsnN5nyiPx4q2c+uatF0hOoJ8hriQwyM9SRIO6k7k6pNTBAaGTGtRU+NBQoZkHrUXduA6wO7AgbCmOa+CAYC7cJuhPKSQ6fw/Xr1xs9PEuCuhWp9utpIZuWGCEUuWnTJkf217Fjx3z++n5nKODvQ5/skMB9XJnN6JFKK7/w6LIoDQZkzR8+fKiYHWQKQl7CDDshxEmwipeijb7sVoie30iBdicDpBeoUsc5FytWzNLV6pYwFBkzZlRWrFgR4v3Lli2jL7IdYUNhnqIs9COQbTAR7F6wYIFpVYvhXoH8NcZqFo0q6CDJhj+FCxemMfqCESNGOGISeB99rYMl+8HrWbtiNXQxFJ06daIVQXDVjviw4T6zBMq0hg2Fb9m2bZvb1e7hw4epMEvuLtCv3azy2uhBDlUDM1UJ41rByOLaNWrUyCeGFq8xYMAAykYygvHjxzukY+ykPmw6QwHdfMhHY0WAAih0lMLx008/0W24D4+xI2wofAdSYfFZwq7BXRYd0jyxSpUpkCgqmzFjhml3F2Zj+/btju5wo0eP1uU10CLXLK4efK6kFDuyzxhFv8ps6LAg9xlfYmSh4MDvuM2sKzotYEPhWzdBmTJlqAmPJ8HPixcvKmXLlnXsLpB550nPaT1BPAJSEmYHvbZloFfrHuSnTp2itNQuXbqYxnhDkRZzFs7Zrqn8puqZjS8B2i5C88kKX4jwwobC93jzucKqFatiqZqKuACC30bRpk0b2hWh9sgqMh+QJ9dS7gLNmPC86FFjhkC+pG3bto7aHLMUH9rWUPgbbCisAdIf1cq03bt3162YKyQwKcrWnHDvmB1cH9n3HNdOy8K0P//803Td9jAe6NXhfEeNGqX4M8/YUGgLGwrfcObMmXBP7CjU69Chg8NYIMMHlcm+BGNAFqBVQOpqqlSpHB3lZEqyt6A2w9cZTWFh7ty5jriWUQWQZoANhcawodAfTE7p0qWjbBwtqnUhrS1XjqjvgT/eLL5yM4JMMlwnXK++ffuGqfIbaq3IQPP1Ls5b8DkoXbo0nStarvorz7yY1yKGW4KQYTRSVv306RMdOXPmDPfz1ahRQ5w6dUoEBASIf//9V3Tq1ElUr15dBAYGCj1Yv369WL58ubAqhQsXFrNmzaLfR4wYIVavXu3V/4cqbZQoUSyhrgwl2SlTppBK8YoVK8S6deuMHpL58Ynpsji8o/BdCqPWbiIEuqFDJtNo0YdZ6/ajyLJClpYdlEqRpYTzQFtjbzMZ0Q7ZSkq/Upgwffr05C70N56x60lb2FBYH0xiMsgMF8m0adM0c0XB1QJ3DYoAze52CQ2kJBcvXtxRuR3S+eDaQSTUyj5+VKmjBgzn2r9/f8XfeMaGQlvYUOgLqrB9ET/A+weftAx0o/BKSwkLI0Tv9ACTv+w6CJ2q4EChLe5HO1tfyYDopcUlhRJ9Lf9uNGwoNIYNhX5gyw83R8GCBanYTm9gkKAsIPWiUK175cqVMD0XMnzsGiBft26dw6AuXrw4yP03b96k5AO0srUyeP9ks6zPP//ctu9ncLCh0Bg2FPqxY8cOih9AdNKXcg8ohkuSJIlD/wcTozdgrMjwMUMluF7AnYbrg74gwa22rZAK6wlYKMiMr0WLFin+wjM2FNrChkJfMNFqHWD2BKyKIT0tGyMNHDjQY2OFFF64K2LEiEG9WOwIAtNSHgUdBytWrGhIHwtfMGzYML/rhveMDYW2sKGwLwjWduzY0eFmgW6Zp9XEqAT/448/FDsDYT/0rpDXB60E7Ch9gc8B4i1maherN1xHwViCFy9eGD0EES1aNDF58mQxf/58ET16dPHXX3+JYsWKiYsXL4b6fwsUKCDq168v7EyKFCnE77//TrUHoHnz5lQvYTfwOUBtBcDn4fjx40YPyVSwoWAM4cGDBzQJNWzYkArijKZp06Zi//79Il26dOLy5cuiePHiYseOHUEeh8nk1q1bwh+4f/8+/SxfvrwYMmQI/f7TTz9RIaMd+fzzz+nziKLPLl26wNti9JDMg+77GxvAriftmTVrFl1T5OybCfRUkXELxCDmzJkTJBMIwW8rtGENrwQKtJDQQRAgdhMQEEDnj3oUMynCah23QtwJ5wlRQzvzjF1PjNlp1aqVOHbsmBg7dqwwE8mSJRPbt28nl9L79+9Fy5YtRb9+/WiVmTVrVlG0aFHafSRJkkTYmV27donXr1+TLAlW1hEjRhQLFy4UqVOnJrdcmzZtbLniTpMmjejduzf9/t1334k3b94YPSRz4BPTZXF4R+F/YAX9/fffO4K49erVo5oPBHKtXn3t6fnPnDkzSOAaWU+RI0emazJ58mTFjuB9lhXbQ4cOVewKZz1pDBsK7UBBk5Um2vnz5zsmxiJFilAWkF1TYefNm+dRwdmYMWMcUihQnbUj0OySUuS3b99W7Ai7nhjTsmXLFpE2bVoxZswYYQWgZItsqJgxY4rDhw9TRtTff/8t7AQWjHC1IaOpf//+oT6+W7duonbt2uLdu3fiq6++Ek+ePBF2o0GDBqJUqVLkfuvTp4/wd9hQMD5lwYIF4uHDh+LmzZvCCsyePVu8fPmSJo0sWbJQxhN+h6y4XUDqK2TZY8SIIYoUKeLR4+fOnSsyZMggrl+/TgbGbvEKnOMvv/xCPxcuXCgOHDgg/Bqf7HEsDruetAM+b8gkeCthbaSvHnpG6OGNQrzy5cvTZyFSpEjkqrEyrm4mb10sqE6H+wnXw+qaTyHRsmVLOr+iRYv6VGLGF3CMQmPYUDBqhdimTZs6gtzw11uRTZs2UbpreJVf0esD1wF6Xb5uOeurNrFx4sShc0S8yk6wodAYNhTa7CSsosx59OhRZciQISH2jsbKsnv37g5jAfE8q5ybFPOTjZYGDx4crufCeX/xxRf0XHny5FHevHmj2I2RI0fS+aVIkYJ6WNgFNhQaw4Yi/Pzwww+0fd+yZYtiZrDCzpQpU6iTKCbIESNGOIxFmzZtQjQsZmTnzp1KixYtNMlAQ5GiVOLt2rWrYjf+/fdfx2eiX79+il1gQ6ExbCjCByZQtCC1SrUrXAxQS/VEHHDGjBmO3hZ169Y1beovDOClS5d0e/61a9c6jOZff/2l2I2VK1c6XGxWia+FBhsKjWFDoWiy6sQW3io9lb0Z59KlSx1B3QoVKphO3gJ+dtSAoIgMv+sFuuHhGkBt1m4SJ58+faL3Vi4I7ADXUTCmA9IYvXr1EpEjRxZm5Ny5c5QzL/FmnHXr1iXV2dixY4utW7eSuFxgYKAwkzLqs2fPxKtXr8SNGzd0e53Ro0dT3QkEHyF9YqeU2QgRIohx48aRlMmyZcuCFYy0NT4xXRaHdxT27iONVXbKlCmpLSr6d4cVVCknTpzY0UsavafNwsWLF3V1PUlOnjzp2F1NmTJFsRsd/9e7JF++fJaKSQUH7ygYU/Dx40dRqFAh0axZM1plmhUU0UEA8O3btyJevHhhfp7ChQuL3bt3k7DchQsXqDDPk74WWoOV/PDhw8W6desct6FYMHPmzLq/dt68ecXIkSPp9+7du4uzZ88KOzF48GARP358cfLkSSrG9Bt8YrosDu8owsb27dsdstxmv3aQl9ZqxY3nypEjh6O1pq9bpc6ePZteG/n/RmhTIX24cuXKjpW3WQP84a0dSZw4sfLkyRPFqtgmmI0exjKTQh7Y0kuQs42tYMKECZVYsWIpX375JQVN1WD7X7VqVdKYRwrfd99953VAlQ1F2Dl06JBp24Xq6TpAMBd1BTK4e/bsWcWXNSuVKlVSpk2bphgFDJR0w6HmxE68e/fOsRCw8rnZylDkypWLfMjyCAwMdNzfvn17yuTYunUryQmgCU7JkiWdJoLcuXNTU/jjx48r69evpw8vCqS8gQ2F/UBRXZYsWZSDBw/q9hr4rGJFjc8OCtz0rFxGDEJd9GcGuYlVq1Y5FnioBLcTGzZsoPOCsvD58+cVK2IrQ4EvWnA8ffqUOpAtWbLEcRu2+Djx/fv3098wDMhxV+8ypk6dqsSNG9erICsbCu/AhBVeaQi9kd3asAvVk0ePHikFChSg18KO9u+//9b8NdA3At+FCRMmKGajQ4cOjqpm9SLPDlSrVo3ODT+tiK2C2ZcuXRIpU6YUGTNmFI0aNXKojh49epQCkBUrVnQ8Nnv27CRhjd7HAD/z5MlDqZmSgIAA8fz5c3HmzJkQXxNBTTxGfTCegxTR9OnTiwkTJgiz8scff4hOnTrpHpBMlCgRSasXLFiQUmaROqu1TPmLFy/ou7Bv3z7TpaT+/PPPIkeOHOLevXuidevWphtfeBg7diylUSNpYMOGDcLOmDOp/X9A+3/evHkiW7Zs9EFDxkHp0qXF6dOnqfF71KhRKQNBDYyCbAqPn2ojIe+X94XEiBEj6LXCy8GDB2msyMFG/rX6iBQpEuXdI8smbty49BNHnDhx6D4rM2fOHJISv3r1qjAruNYTJ070yWslTJiQjEXlypXFkSNHyFjg7/z582vy/F27dqWFVM2aNemzZibQx2PRokXUQnbVqlVixowZol27dsIOZM2aVXz77bdkMPAenDp1SkSJEkXYEsVCIMMAbqNZs2Ypv/32G+Vru4IK1F69etHv0N9B9oVrm0OcNtxSIYEsDWzH5IHc+rC4nnr06BEkGO/JgTgKXG4Iwrdt25Y0h3DOiMXcuXPH9AJ0CPZBBiM8NQl6gGuIKmojP7/Qu8J7jAQMxEnCwr59+5SGDRtapsod/Pzzz3TeSCrxdRaYnmzevJkk53Fu48ePV+zqejL1jsIV7B5gxS9fviwqVapEHbaePn3qtKtAvn7y5Mnpd/w8dOiQ03PIfH75mJAqWXGEF1SpoiHMp0+fghwfPnyghjiomIVrCz/h8gKPHj2iA7nawYFdB3ZZcLXhZ65cuSiHH43vzbCixKqqadOmwkygOx1Wsqjt2Lt3ryhZsqTPx4DP6aZNm8j9id0m3KabN2+mWhNPwWcGn6nHjx9TzYJVuq+hKx7cM9hJff311+QW1uI7ZgSKoji+Z/juSQYNGkTu8cSJEwvboVgISPwiJx95zDKYrV4hIvsguGD2gwcPHI+ZPn067Uq8ye32VTAbY8JYUd2KsSNIOWjQIKV169Yk5YwsHbl6Ce5AGiYCa0gCWLNmTZBUYb1BSqhZdztYfXfu3Flp0qSJ4WPE56hEiRL0nsWPH59SiL1h+fLlyldffWU5yWvshhMlSkTn3b9/f8VqnDlzRqlfvz55KtRABBHZlTgvBO+tgm2ynuC62bFjh3Lt2jVl7969lOYKt4wUHEN6bNq0aZVt27ZReiy+fDhc02Phfjpx4gSltCHzxMrpscjWQk4+Josff/yRmujATRWSAcmePTt9eFHLoDaYWoPJt1ChQuRa8WXNgLeYRXYBwoGlSpVyFCS6y4aCq8Yu7hpkKcoOgd4aSKPZt2+fQ0EWHQ/VYJ7CfViYWqWBk20MRYMGDSitDrGIVKlS0d+XL18OUnCHL1rMmDGVOnXqBFHHvH79ulKlShXyjcLIwPjYseAO6ajYSaElZbNmzaj+JEKECEEMB26Hyue6des0bTKD1Rbeg9ixY5tGORQ7NMRKjN5BuDMWcmeB3SBqIVxB/w4Uk6LAy2o7iJDA9xjnDCl3szY6wqJq6NChQYoWhwwZQjv+4ECqNc6revXqihWwjaEwC1YwFMGBfgrQ0e/SpYuSN2/eIEYDE3utWrUoyKuFmwpGGgbIDMA4NG7cmM4T529WsDKVRXnYHUP+w3XCgmAhenWbxQCHF9RTyA57ffr0UcwI+qFjfOij8t7DheWFCxccO3vsMMwOGwqNsaqhCK74a9myZeSyk42E5IHdR7FixahnhJlUT8MDdldwEyAzxczAGECaBu9D1qxZlStXrjjdj120lTKcPAGuU+mqOXDggOFZenDNqid37HTgiVi4cKFX114WGMIFa9adrIQNhcbYxVCowYf42LFjlHpbuHDhILuN0qVLk0x0aNW0cHkF5zIxUwDVCmAnkS5dOof/Hh3j7E6jRo0ccTQjK/mHDRtG4yhbtmy4n+v+/fvkKrRCN0c2FBpjR0MR3IQKeRN8WdSxDUxaWFmhbiU4fzLy4/EYZFqZAQRIzer3Dg2o1yKWJgPcdolJuHONQl0X5wuxTl8tkBCUVisFo94HMdABAwZokuwwaNAgOif02TZzPxY2FBrjD4ZCDb44MAAFCxZ02mWgSKxbt25OWU3IusJ9kLY2GkwAiLuUK1fOdO1IPQUihdJYfP7555Y1ep6yevVqh+sTmY1607t3b3o91xRXLbPhXrx4QckJeB24P80KGwqN8TdDoQa1KVhpQaXX1TX166+/OrKtzOBDR/8L9GBAzYmZV3Lq1e2MGTOCSFVjV4TsMVznGjVqkA/dziBLD+eKOiEoJ2gJ4j3qnhE7d+6kuBWyJfVk6tSpDpUFs84bbCg0xp8NhXrFhYwmZEmpazawy0BdilliAadOnVJevnypWAFI38vr6Jolg7+jR49O90Guwyz1H3qAiRyZXTjXrl27ava8MMDYqYwePdrJOMPlpTfv3r2jxASc0/fff6+YETYUGsOGwpnbt29TLrw6cwpV8nBDobDRlyDY7lo7YyVQoYyJLLj+ETDM6HcgXSVmz6IJD1AikC4orPrD2ixJvbOFskFwbiZfZ3bFiBHDNAspNWwoNIYNhTPwJcu8/0WLFpEbSu2WqlChAk1yejfPQRwCGVsZM2Z0KsQ0K1hlIv3YtarXHcicQQoprusPP/yg2JlWrVrReeL99HZXiJRvGFXUDUngxjLyc/Hp0ydqpGaksXIHGwqNYUPhDFRs06dPT19sdRD2P//5j5NbCkV+kGzQy2BA2iVDhgykH2QFiQtcH1wX6DR5sztAHENeU2iV2RXot8lYWKdOnUJ9rJqePXuasjXpnj17HPUiZpO2YUOhMWwogoIgdnCN5SGZApkUBJXl5Aa9LRQ06WEw4G4Iq1y3r0GQGhXJixcv9vr/IqFATjgQfLQraJkqPzfQcHMFBhbCjpD1UWsqwR0akrSG0dSuXZvOp2bNmoqZYEOhMWwovAcBQ7hKoNQrv/jQ9vn999/DFZiF+8YqomvIuMHuS01YC8swQbZs2dIhvaJnr2+jadeuHZ0ndq1wL7pmQsmJF248K3Du3DnHTnvXrl2KWWBDoTFsKP4L1Gq91XKCPx7FePHixXMyGFgVexuchZGoW7cu7VZ8kXMf3t0DKnThFtMq2I7zR+qvTLtUF43ZCRgHWaWOojVIsatdTdg5QC3aisavePHipklKYEOhMWwoFMpmkkFVpHV6C9xUkAvBl14ajDJlynil84OVJYrQ4HYwi/hgSKCOo0CBAlTp7ir0F95iLsi5y0lUT+l4I5A7LuzE1AkSSJqwMnfv3qWdIM7FyC6LathQaAwbiv9eA8QeoM8THrDDQKtaFD3JSaBevXoe60XBWOzevVsxG1glospYvVrEJK5HXAZ6Qgji49qh9a9V6kbcASHKgIAAMn7ymkEOH+eIKmc7fPcG/C/OhMJCMxRRsqHQGDYU/49W22assps3b+7QlUJqIyYGKNyqwRcKPRnMfk2gh4XzgDy1ryrmZbc4dDU0Q2W8t6hjD9hJSPfk4cOHHbsnpMriNrhu7OBSS/o/efXJkycbPRw2FFrjz4ZC75UPOrtVrVrVqdIb8gcIeOO10YwKxgQNiMwMAqtwLbg2utFb20pWb6Ndrll8357Eb1D/go6ValatWhUk7gJZFvnZMONO0lsmTZpE5wKDYbQeGRsKjfFnQ4FaCeT9hyY3Hl6QCpknTx7HpAD/PnLQoe8PNxX6EpsJrHrVfSOwokdqsK9ZsWKFI3aEGJAZgcFXp1Kj/gXjRazJk8+VLMSDHLk3ve7Nei2yZMlC5wNXlJGwodAYfzUUqGqFSyg8sgregMkWapvqgDe61AWXT28kkIbA5IzAuhlW8XBjyOs1d+5cxUwgCI3VMwy+GtTVeNqxD3EtqcZqFjl7LfqGIysOQW6jYEOhMf5qKADSEMeOHeuz18PkARE1rCJl/ALpsOPHjzeNMB4MKPR7INantdppWEFLUam55QujHhKY+NS7B3QXxLjQwS88RhWGRZ6f2SqcvQXXAd0kjY69sKHQGH82FL5OKc2RIwdd6+HDh5MvGy0l5WoZvyOm4WvgYnKtpob7xEwgUwguQlwnBLld26n6AmTFYaelXlhI1eHwxrowuSJoj/MrVaqU7jpierNr1y46FxTiITHBCNhQaIy/GQoUsxmldgktIxRbSe0mTAgIEMsKb7jCoLjqq4Y+Z86coRgJjqtXrypmBrsbWWOBokZXPSSt3YSIG6kNgHSBqTXAtE6hlW1GfZk0oBeQ9MC5IGHDCNhQaIw/GQqkp8IfjFacvqp+dXVJBNcCFFo+UroBB7T+feFiwdgQi0CGjhUUanGdZG8HVHHrkTaLa5IvXz56DXVvbxgmva8RXJB4XSwczCjd7e0iRCYiIHHD17Ch0Bh/MhT4oqMFKlakvli1I9iJSdjT11q2bJmjz7L08WrZWxo7GUhCqzvk4X03Q9Dam4ws2U41vI2AcN6oyh83bpzT7XheyIj4qm5E7cqS7sgvv/xSsTqtW7emc4Ecua8/Y2woNMafDAWAO0FL2Ql3woEywwkrRU9BsBSTuTQWKMrSIsce5y2bMam7olkR9LGQ1wcy5eF5j6Sg3YULF5zeA6Oqi6H1JBs6IT3Yyty5c8dh1NHoyJewodAYfzMUvgSprwiChiU4if+L5kl4b5AhBWkQb/PsXV0zSC9Fn2oruJlCY8iQIY64DgrXQgPihT/99FOQBknwoWP1bibVXpnllSpVKst/L/v37+9wp/rS+LKh0Bi7GwpM0qhXwMpT7+0vVqhapjfCLw4pEHXvC09FC5HPDm0hdazDSi6m0MC5yGZJqHgPTm1WnXIMl5WUMVfrR5nxmkDyA+8dxtuxY0fFyjx79ozceDiXWbNm+fR12VBoiN0NhXRToFJW7V7QGvj/IWaHrCZXTafwAhdEkiRJHLn2kEQPre4C7TPxeOwg7AomVOnTR2WzrHHAe44OhOrdAwzC119/TRlFZqkPcQc0wORu0uyy86GBlGKcC3bIvqo+Z0OhMXY3FJggULcwZ84cXV8HFbaIJ+DQo3Up1Fpr1arl2F1A4htZQBIEZaG8qi4OGzp0qKbBcLP6waUYXeXKlcnd9uuvv9LfkE2xMnI3mStXLqcEBCsa9JT/y1aDOoEvYEOhMXY3FHriGntAcREMhp5GD3GG2LFjO4rPkMIJQ4iVp6uUhD8gmx1hx4if3377LbnssDCAK9DKYGcqd5KIyViZKVOm0Hkgq88XOzpv5rWIgvFLnjx5IkaPHi0+fvyo22vcuXNHlC5dWqxfv95xW7Zs2USCBAl0e80IESKI5s2bi6NHj4oCBQqIx48fi+rVq4sTJ05gUSSeP39OP+3Iu3fvxMqVK0WXLl2czrFQoUIiRowYokGDBvT3hAkTxNKlS0WLFi1EwoQJhZVJlCiRGD9+PP0+bNgwceHCBWFVWrVqJdKnTy/u378vpkyZIkyF7mbLBthtR4FVt1xltm3bVrfXQRYSXgNxCV9mcyBjCf2lsYXH6lkd6LZDNpMatbsF7gvZRe3o0aOO27GDkytUuNpwPyrNIZFit88zuiZaWd5j7ty5jp2w3jLk7HrSGLsZClnohg9jWNqaegqCcpBz8PXkLKUk0qRJQwYKfQ6Q9SMFBl11m6wIelGgMNK1p0Pnzp2pGC6ka45JVMZxcH08VXA1O9DekkYS6r5W5f3795Qm6wtXGhsKjbGjoQBaB3FlHr6v0ylRgIVAtQRV3kj3xWQqQQHhZ5995thdYDI1QztKT8DKEhXpakFEZKfJnYG3/mzEJ+RkBHkSK3bHC44xY8bQOaGIE59Fq/L77787ZEr0jCGxodAYuxgKTC569VfGZJUiRQq6Tsio8RVI5ZQTXmhgQuzbt6/DWJQuXdqUE4qroZVV6HCjqcHOSJ3F5a3OkBTY++677xQ7gPdXiiLWr19fsSofP350NPHq16+fbq/DhkJj7GAokDoKVwPE3PSS5xg2bBh1ptPTnYUvkdp3C0VRrKobNGjgsV4Uai6kGi2MmxGCbCEZiBYtWlDWC85LPV50RUO7VT0a6OBAvwc7cOzYMYfQ3qZNmxSrsnLlSjoHGHN8d/WADYXG2MFQoD8BJiC4HNSNZcIzqaGSW93+E64cPd05SHNFNS4K5dSExc+O3RVy76XExYQJE3zqMkMQGrIav/32m9PtJUqUCKLPpOe4ZMIBJqRTp04pdkAmMOCzbtXWqZ8+faK+4jiPbt266fIabCg0xg6GAsDHrZV8BvSZZAGXnhOZ+rmh7YTXhHCfFgYJMRrsROSqulGjRrrlr8Mtoo4JHThwgF4Tcu7qCnJUG+M8fTXBYVwIiGMs2LVosYgwGsRgZOtUVOhblQ0bNjjiUOrCUa1gQ6ExVjUU0PbRyw2EFTkyiSA/rUc6IlbbFSpUcJK3htFYuHChpnEWPCdeQyqkwset9Zdy1KhR5OoaNGiQ0wSdOXNmCrobXfQWGBjoEFeEnImV00slsvIcyqxqN56V+PTpE8XRcB56FIqyodAYKxoKSGRgVYWqVS1aLWInAp+2Gr0C4+ogNdwHvnAJ7dixg9KFZdwiLDUGGCdSc+vWretUfS7PpUqVKkEebxbQpAorV4xz8ODBitVRT7JW7luxc+dOh36Z1h0W2VBojBUNBcaKPHsEr8OaGaMOEEL+AXnqevRiRjAZmj3YaquNEDI+1DEQvcEXUcYtokeP7rbeAoFz9D1GjYYaCO/h/yMYqY6hYCIOTaTQLMVekDpRd66zsqtV7hTRttWqVKpUic4B3xEtYUOhMVY0FHKCCqvPWb3ahSuiXLlyStWqVTVxy+C51c8vU1YDAgIUo8F7XK1aNUfcAuqqOH+4Z9S7hI0bNzqqztXAjYXsr+Akva0AXByyFsGq56AGgWCcD9x8vuqzrjUynoVsLi28AxI2FH5qKJYuXUqFWeEBwVzIPMBXrw4Ya9UOFP76bNmyOcVOkNOPlqZmSVNFRhJqC6SxgDHAz19++cXpeqBpTr169SybWRPSucvMK8iQQxbEyuB9kvU9+FxblRo1atA5oL+IVrCh8ENDgVagWHHAl6nW+QmLoZAZI5D50Lr6G/57PHfPnj0VM2bLoHcxAs9YfcIVg+spDQb6G/sDallyPbXAfAU+x9KdCKkPK3LixAnH5xBKBFrAhsIPDQX836hGxWTmjS8cGTeYEF0zRvDlCo9PHTEGBG8R11C7bLBrmDdvnuHXEj5rCMmhDaUEOyYpWb1//37HeKVOVPr06XXpo2FGUKyGWAXO27XWw2rgfYXrFOcCnSurUr9+fU3PgQ2FnxiKW7duOcUM4Cryxj2E1b6cBOEHDSt4Tfizt27d6nS7lCFwzZbyNcjiQfxD3b1P6umg+5sanAP8wOoUUQjswcct6x4QxPYHEJ/BOaO3h5a+cSNAv28UVuJ8rBqoP3v2rKPqXAvlXzYUfmAoIOeAD80333zj1f9z9adDMgK+aLiuvEFtkGQKH3zB6tvxnL6cYBDrgLy4aw9l6XNXu9LQ3Q6prJ666ZAYULx4cUcBFFqJ2h3sKOVK3A7xChl3QodFq55L06ZNNUv8YENhQ0OBFFe1xtHmzZtpTDVr1vRoF4HHDBgwgFwrau0Y7Cq8KbCC2whpt1DqlMCfHy9ePGo9qmf3OjXTp0+nAB8qmSWHDx+ma4JG9eprAgOBWobw+qcRv6ldu7bDV/zzzz+bqhZCD2BQZbwCCQdWBt8fJCDgXNTFj1biypUrjp1ReHe2bChsZihQT4CgqjrrBpN7SKt1jBPZT2itqEZKNaCRe2hgAkRm0ujRo51WX7IJPKQ71GjVrxirWHUaIyZ3GCAYJzVNmjQJotmPccJdguwvvaqLMT70fJDGAr+bvT5Cy3iFFgkORgLxQ7kr1KMmyBfAYMsmTeFZqLChsLihQIBZPfnIRjyYHIMLGsP1ozYacMHg8Qgkq3sNQP5i3bp1IX641BlKeIxcfaFmQIICuAULFtBKM6zgubHzQCGfGmRC4QusVknF4+SkrN5RrV+/ngynEUJ2GD92E3JcCDJqZSjNyvfff++IV6hjPVYD7x2kYXAu1atXV6zIrVu3HFX04VHIZUNhYUPRpUsX+hCoK0nxuljdY5XsmnUjVxfqngJ4HNI8kdroWnCHiRWrKrVhQQUrgtpoFaqmU6dOVGQXni3u3r17qQANBkqdhiorgNUifHAHBJeGigppBO/M1mgI45Lps/AZ6ylpYjRYcGBnh3NFtb9VffwA3yH5vq1evVqxIl26dHEkY4R1V8GGwkKGwvULh85reC1IJatXqfhdZihhRSGZM2cOrfwhF61eceMxEJxD3YIaqZaqjjFgd4DbIHfw6NGjUMeMVT7iAWpjA3cRPrSQMldPmAMHDgw2Hx+1Cgh+q1t2YhyQ0bCSKweyI7IFJ4yzr2I0RtVXyPRhq8cr+vTp40h51ksxWO+YpfzchdXYsaEIgUmTJinp0qWjFTsmtYMHDxpmKLDqh1wCejhfvHjRybWD6mV8gOvUqeP0f7Dih+YS/OLq50HmkmvsQU7+yIxSu59GjBhBE9rEiRMdfnz8xIcN9RTq9qEAY2zYsKHHkz9uV+96sDPC+OD2UmNVOYXgwDWD5AXOHSnB4XHLmR24IWW8AinGVgWfZzTykjItVjZ2efPmDVNMjg1FCG4CTLJYgcOHj/aS+HJ70j1KK0PhulKGWwfPi1W4un2oTDfFeGfPnu24Hamc0mWj7ishe0MUK1bMkf6KDw7kJdDoB64fCbap8ouubgMKAxKc8Fhwkz8MSsqUKWn3owaZWBijnYyAp8B9h/dRpl9aNVDqbbxCvcixGkh6kN8zK+paPX782PH9DEuHQjYUwYAdhLrmABMpJjtMkHobCmxtW7VqRdr40A2SvY/hd8fkLoOiSD2V7ijsJqQ7SE68GDN2RLhdXSsAd5F8DrVSLJq24DbUFqjBLga3q78cqAtAzrzr9cAuBYedV8laAeMAIyFrSuzSMc4V7FCRcSPjFVZdGGDRhOw9nAdUBKyY6jzof3E96KepPQeewIbCBfj3MeGi97Br8QrqEFzBqhwXTx7w94fHUOADKIXl5CFX/qgDkCt8tSsHRkQ+Vt3qE5k+mIxc88DhIsIKX90EB9kpkMF2nbC8rZ1gPAcGFS5CWcWtjifZNV6hR1MdX4EdEXYUOA/X+cEKPHv2zBG7nD9/vtf/lw2Fy4caF8TV/450TFcJB7UP3vUIj+sJMQD0wMVWETUB6lRUyEYgZoAMITVwkaE7l5WCu8x/g/2oBofhtuIq1Zt4BXbl6EZoZaD3hfcLwntWZOTIkTS3oPmWXoYiAv4RNufu3bsiVapUYt++faJEiRKO23v16iV27twpDh486PT4t2/f0iF5/vy5SJMmjXj27JmIGzeuT8fOWJM3b96IaNGiiYgRIwq7n2eMGDGElXn37p2IHDmyZd+r9+/f0/gjRIjg1f/DvBYvXjyP5rXIwg9InDixiBQpknjw4IHT7fg7efLkQR6PLzgOhgkrVp88/ek8o0aNKqxMlChRdH8Na5rQMHwQChUqJLZu3eq47dOnT/S3eofBMAzD+OmOAnTv3l00a9ZMFC5cWBQtWlSMHz9evHr1SrRo0cLooTEMw5gavzEUDRo0EIGBgWLAgAHi/v37In/+/GLDhg0iWbJkRg+NYRjG1PhFMDu8eBP0YRiGsdu85hcxCoZhGCbssKFgGIZh3MKGgmEYhnELGwqGYRjGLWwoGIZhGLewoWAYhmHc4jd1FOFBZhAjnYxhGMYOyPnMkwoJNhQe8OLFC/oJYUCGYRi7zW+op3AHF9x5AHShoEAbJ04crxQapersrVu3uFAvBPgauYevj3v4+oT9GmHqh5FImTJlqMq5vKPwAFzE1KlTh/n/483hD7F7+Bq5h6+Pe/j6hO0ahbaTkHAwm2EYhnELGwqGYRjGLWwodATNjwYOHMhNkNzA18g9fH3cw9fHN9eIg9kMwzCMW3hHwTAMw7iFDQXDMAzjFjYUDMMwjFvYUDAMwzBuYUOhI5MnTxbp06cX0aNHF8WKFROHDh0yekimYMSIEaJIkSJU6Z40aVJRu3ZtceHCBaOHZVp++uknUgTo2rWr0UMxFXfu3BGNGzcWiRIlEjFixBB58uQRR44cMXpYpuDjx4/ihx9+EBkyZKBrkylTJjF06FCPdJ2Cgw2FTvzxxx+ie/fulJZ27NgxkS9fPhEQECAePnwo/J2dO3eKb775Rhw4cEBs3rxZvH//XlSuXFm8evXK6KGZjsOHD4vp06eLvHnzGj0UU/HkyRNRqlQpESVKFPHXX3+Js2fPijFjxogECRIYPTRTMHLkSDF16lQxadIkce7cOfp71KhRYuLEiWF6Pk6P1QnsILBqxhsl9aKgt9K5c2fRp08fo4dnKgIDA2lnAQNSpkwZo4djGl6+fCkKFiwopkyZIoYNGyby588vxo8fb/SwTAG+Q3v37hW7d+82eiimpHr16iJZsmRi9uzZjtvq1q1Lu4uFCxd6/Xy8o9CBd+/eiaNHj4qKFSs66UXh7/379xs6NjPy7Nkz+pkwYUKjh2IqsOuqVq2a0+eI+S+rV68WhQsXFl999RUtMgoUKCBmzpxp9LBMQ8mSJcXWrVvFxYsX6e+TJ0+KPXv2iCpVqoTp+VgUUAcePXpEPkJYdDX4+/z584aNy4xgpwXfO9wIuXPnNno4pmHx4sXksoTriQnK1atXybUC926/fv3oOn377bciatSoolmzZsLf6dOnD6nGZs+eXUSKFInmo+HDh4tGjRqF6fnYUDCGr5pPnz5Nqx3mv0AOukuXLhS/QSIEE/wCAzuKH3/8kf7GjgKfo2nTprGhEEL8+eef4rfffhOLFi0SuXLlEidOnKAFGSTFw3J92FDoQOLEicmKP3jwwOl2/J08eXLDxmU2OnXqJNauXSt27doVLhl3uwG3JZIeEJ+QYEWI64SY19u3b+nz5c+kSJFC5MyZ0+m2HDlyiGXLlhk2JjPRs2dP2lX85z//ob+REXbjxg3KOAyLoeAYhQ5g+1uoUCHyEapXQPi7RIkSwt9B/gSMxIoVK8S2bdsohY/5fypUqCBOnTpFq0B5YPUMtwF+93cjAeCqdE2phj8+Xbp0ho3JTLx+/TpIMyJ8bjAPhQXeUegEfKew3PiCFy1alLJVkP7ZokUL4e/A3YQt8apVq6iW4v79+44mKsjK8HdwTVzjNbFixaJ6AY7j/Jdu3bpRwBaup/r161ON0owZM+hghKhRowbFJNKmTUuup+PHj4uxY8eKli1bhu0JkR7L6MPEiROVtGnTKlGjRlWKFi2qHDhwwOghmQJ87II75s6da/TQTEvZsmWVLl26GD0MU7FmzRold+7cSrRo0ZTs2bMrM2bMMHpIpuH58+f0ecH8Ez16dCVjxoxK//79lbdv34bp+biOgmEYhnELxygYhmEYt7ChYBiGYdzChoJhGIZxCxsKhmEYxi1sKBiGYRi3sKFgGIZh3MKGgmEYhnELGwqGYRjGLWwoGMZkNGnSxKGKaiTFixdnkT2GYEPB+D3Nmzenvt2u7Nixg3pVP3361PF3rVq1SLkU2kvoOAcpZ1f++ecfknSGQB0EIiHtDI2dmzdvhjoWNJhZv3499VYIC4MGDaJxBac9hsZQ6LLoOuYlS5aQNpAr33//PSmQhlVIjrEPbCgYxkP27dtHvauxyv77779J4LFp06Ykla42EliJb9myhXojXL58mZoQ4Sda46LhjjvQ0xhd22LHjq3ZuNesWUMijJs2baK+ya1bt6bmWrK7YP/+/cXkyZOD/D90Q3vx4gX1pGb8HK3FqBjGajRr1kypVatWkNu3b99OYoVPnjwJ8f9WrVpVadGihePv9u3bK7FixVLu3bvn9LjXr18rqVKlUr744osQn+vDhw9KvHjxlLVr1zoJS+bKlcvx94oVK2hMU6dOddxWoUIFEnyDqGJwQosjR45UGjRo4Hh80qRJlUOHDtHvbdu2VcaOHRvimHBujRs3DvF+xj/gHQXDhAOsyGWvb7hosHtA3wjXBlWQT+/YsaPYuHEj7TqCA7sUPB+k6SVly5YVZ8+eFYGBgfT3zp07qTEW3GDg/fv31Ie9XLlyokGDBqJHjx4kK33v3j06cFu+fPnEkSNHxJMnT6gp0ps3b0TmzJmpqyDarbpzc0Eif/fu3ZpcK8a6sKFgGCHIfQR3j/oIrRE92k2iV7PsMYLJHPEMdFoLDtwOsWa4oYIDHcjQXCZp0qSO29B/AoYIBgLAQMAYyL/RhwHGAr0ZYIww7siRI5OhwoHbAgICROPGjcn1hXjM/PnzKcbSoUMHco+h93S2bNmoGdCZM2ecxoT4ClqzcpzCv2FDwTBCiPLlyzt1lMMxa9asEB+/fft2MhAzZ86kFbyasCr3Y6UfLVo0CqBL8HuZMmXIQMAIYXeBnQnaoZ4/f54MBgxAzJgxQw1yw0Chc16dOnWoJWbFihVFlChRxLBhw2h3gdgFYi5qYGhgJPB6jP/CHe4Y5n8d5OCOUXP79u1gH4vJGVlC48aNc5pYkyRJIuLHjy/OnTsX7P/D7Zj4XV9HApcSWli+e/eOsqUkcCuhcxtcQAUKFBBx48Z1GA+MBe4pb4CBWbhwIXU9mzNnDj0Xxo5OccjOQgAbXfYA3GS4Ntx50L/hHQXDeAEm52rVqomRI0eKtm3bOt2HHsWYbJFhJNu7qncLU6ZMITeQjGm4ItNasWtQI+MUSGOF0QD4icyqvXv3Om4DMDAfP34McfzY7bRr147aYsJNhcfCdQXkT/X/P336NBknxr9hQ8EwHgJ3E4wEgr9169YlY4BDHZxGoRxiA5UqVaK0Uvj3d+3aRQYCE3FwaagSrOoLFixIbiA1SMlNkCABGSC1oVi5ciW5hBBbkKRPn15cu3aNXGdIgXV1GcGdhteRdRP4v9u2bRMHDhygHVLOnDlpVyTBLqZy5coaXD3G0hiddsUwVkmPxeOC6/WNftZqAgMDlc6dOytp0qRRokSJoiRLlkxp3ry5cuPGjVDHMmXKFKV48eJBbsf4IkeOrLx48YL+/vjxo5IgQYIgj/3333+VunXrKvHjxw/Sh/z+/ftKunTplDt37jj9n8GDBysJEyakvtMHDx503H779m0a/61bt0IdN2NvuGc2w5gIuKiQgfTHH3+IEiVKGDqW3r17U0ot4iOMf8PBbIYxEQgaL1iwwFE5bSRI04X0B8PwjoJhGIZxCwezGYZhGLewoWAYhmHcwoaCYRiGcQsbCoZhGMYtbCgYhmEYt7ChYBiGYdzChoJhGIZxCxsKhmEYxi1sKBiGYRjhjv8D689PF4i20ekAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 400x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, (ax1) = plt.subplots(1, 1, figsize=(4,4))\n",
    "\n",
    "# Plotting results\n",
    "for P in [1000,2000,3000,4000]:\n",
    "    df = isobars[isobars[\"P_bar\"]==P]\n",
    "    ax1.plot(df['H2O_wtpc'], df['CO2_ppm'], '-k')\n",
    "for XH2O in [0.2,0.4,0.6,0.8]:\n",
    "    df = isopleths[isopleths[\"XH2O\"]==XH2O]\n",
    "    if XH2O == 0.6:\n",
    "        df = isopleths[isopleths[\"XH2O\"]<0.7]\n",
    "        df = df[df[\"XH2O\"]>0.5]\n",
    "    ax1.plot(df['H2O_wtpc'], df['CO2_ppm'], ':k')\n",
    "\n",
    "ax1.set_ylabel('CO2 (ppmw)')\n",
    "ax1.set_xlabel('H2O (wt%)')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "volfe_dev",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}